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We propose a dynamical model in which a network structure evolves in a self-organized critical 
(SOC) manner and explain a possible origin of the emergence of fractal and small-world networks. 

Our model combines a network growth and its decay by failures of nodes. The decay mechanism 
reflects the instability of large functional networks against cascading overload failures. It is demon¬ 
strated that the dynamical system surely exhibits SOC characteristics, such as power-law forms of 
the avalanche size distribution, the cluster size distribution, and the distribution of the time interval 
between intermittent avalanches. During the network evolution, fractal networks are spontaneously 
generated when networks experience critical cascades of failures that lead to a percolation transition. 

In contrast, networks far from criticality have small-world structures. We also observe the crossover 
behavior from fractal to small-world structure in the network evolution. 


I. INTRODUCTION 

Complex systems consisting of discrete elements and 
their pair interactions can be described by networks. 
Many of large-scale networks representing complexity 
of the real world are known to have common proper¬ 
ties in their topology M, such as the scale-free prop- 
erty [i L degree correlations 0 , 1 , or community struc¬ 
tures 7|- I n particular, structures of real-world net¬ 
works are classified into two types from a viewpoint 
of the relation between the number of nodes and the 
path length, namely small-world structures [i| and frac¬ 
tal structures Q. For a small-world network, the av¬ 
erage path length ( l ) is extremely small comparing to 
the network size N and increases at most logarithmi¬ 
cally with N, i.e., ( l ) oc log N. Numerous real-world 
complex networks possess the small-world property E3- 
HU. On the other hand, a network is called fractal if 
the relation Nb{Ib) oc l^ dB holds, where Nb(Ib) is the 
minimum number of subgraphs of diameter less than Ib 
required to cover the network and <Jb is the fractal di¬ 
mension [g]. Since this relation at Ib ~ (Z) suggests the 
power-law scaling ( l ) oc N * 1 / dB 0 , the fractal nature 
seems to conflict with the small-world property. Never¬ 
theless, real complex networks that are small world in 
the sense of (Z) oc log TV often satisfy the fractal scaling 
Nb(Ib) oc Zg dB , as observed in the world-wide web, ac¬ 
tor networks, protein interaction networks, cellular net¬ 
works 0, power-grid networks p~5| . and software net¬ 
works [16l 17]. This apparent inconsistency can be rec¬ 
onciled by taking into account a structural crossover from 
fractal to small-world scaling associated with the change 
in length scale jl4j . 

It is well understood that the small-world property 
arises from the existence of short-cut edges Q. Only 
a tiny amount of short-cut edges added into a non¬ 
small-world network drastically reduces the average path 


* akitomo0416watanabe@eng.hokudai.ac.jp 
' s.mizutaka@eng.hokudai.ac.jp 

I yakubo@eng.hokudai.ac.jp 


length. In contrast, the microscopic mechanism of the 
emergence of fractality in complex networks still remains 
unclear though fractal networks and their relation to 
the scale-free property have been extensively studied 
dm em. It is thus also not understood why there 
exist small-world and fractal networks in the real world 
and how fractal networks crossover to small-world ones. 
In order to deal with these problems, it is significant to 
remind that many conventional fractal objects embedded 
in the Euclidean space are formed by dynamics exhibiting 
self-organized criticality [25U3~it . In self-organized criti¬ 
cal (SOC) dynamics, a system approaches spontaneously 
a critical point without tuning external parameters and 
fluctuates around the critical state due to the instabil¬ 
ity of the critical or near-critical states. One of the re¬ 
markable features of SOC dynamics is that stationary 
fluctuations around criticality are accompanied by inter¬ 
mittent, avalanche-like bursts of some sort of dynamical 
quantities, in which the avalanche size distribution obeys 
a power law. It is natural to consider that fractal complex 
networks are also formed by SOC dynamics. 


SOC dynamics on static complex networks have been 
extensively studied in previous works [U l32l44l| . How¬ 
ever, for the construction of fractal networks through 
SOC dynamics, we need to consider the interplay be¬ 
tween internal dynamics and the network topology [HE 
[H|, by which the network evolution itself displays SOC 
characteristics. There have been many models of net¬ 
work evolution driven by internal dynamics related to 
self-organized criticality, such as models based on the 
Bak-Sneppen dynamics ]4bl - t47] , other models of ecologi¬ 
cal systems [48j , models related to the sandpile dynamics 
fril - oTI , a model describing the motion of solar flares [12, 
and rewiring models based on state changes of nodes or 
edges |53 h 56| |. Although these models generate nontriv¬ 
ial networks through the couplings to internal dynamics, 
it is difficult to say that fractal networks are formed by 
SOC dynamics in these models because of the lack of in¬ 
termittent, avalanche-like behavior of network structures 
liM 153, 541, the necessity of parameter tuning for criti¬ 
cality 23 55}, or the absence of fractality in generated 
networks Ea sum EH Hi- 
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In this paper, we present a model of fractal networks 
formed by SOC dynamics. Taking into account the evolu¬ 
tion of real networks, the network instability required for 
SOC dynamics is realized by overload failures of nodes. 
In general, when the size of a functional network becomes 
large, the probability that all nodes in the network can 
escape failures decreases. Failure(s) on a single or a few 
nodes can cause a cascade of overload failures, and the 
network decays into smaller ones. Our model combines a 
network growth by introducing new nodes and its decay 
due to the instability of large grown networks against cas¬ 
cading overload failures. It is numerically demonstrated 
that the present dynamical system exhibits self-organized 
criticality and the network evolution generates both frac¬ 
tal and small-world networks. Furthermore, the crossover 
behavior from fractal to small-world structure has been 
observed in SOC dynamics. 

The rest of this paper is organized as follows. In 
Sect. 2, we formulate the model combining a network 
growth with cascading overload failures induced by fluc¬ 
tuating loads. Our numerical results are presented in 
Sect. 3. In this section, we show the time development of 
several measures describing the network structure, SOC 
character of the dynamics, and the fractal and small- 
world properties of networks generated in SOC dynamics. 
Section 4 is devoted to the summary. 


II. MODEL 

A. Network instability — Cascading overload 
failures 

In the present work, the instability required to con¬ 
struct an SOC model is realized by cascading overload 
failures in large networks. Our daily life is supported by 
various functional networks, such as power grids, com¬ 
puter networks, the world-wide web, etc. Functions of 
networks are achieved by some sort of flow which plays, 
at the same time, a role of loads in the network. The 
load on a node usually fluctuates temporally and its in¬ 
stantaneous value exceeding the allowable range causes 
a failure of node. This overload failure may induce a 
cascade of subsequent failures which reduces, sometimes 
greatly, the network size. Such cascades of failures pro¬ 
vide the instability of networks. Recently, the robustness 
of a network against cascading overload failures induced 
by temporally fluctuating loads has been studied ]57f, 
This study employs the random walker model proposed 
by Kishore et al.[5a |59|, in which fluctuating loads are 
described by random walkers moving on a network. Since 
our model is based on the study by Ref. [IV] we briefly 
review this work. 

In the random walker model [H|-[6Cj, we consider Wq 
non-interacting random walkers moving on a connected 
and undirected network with Mq edges. Walkers on a 
node represent the temporally fluctuating load imposed 
on the node. Since the stationary probability pk to find 


a walker on a node of degree k is given by pk = k/2M 0 
0, the probability that there exist w walkers on the 
degree-fc node is presented by 

h k ( w ) = (^^jPki 1 - Pk) W °~ W ■ (1) 

This binomial distribution gives the average load on the 
degree-fc node as ( w)k = WoPk and the standard devi¬ 
ation (7k = y/{w)k{ 1 — Pk)- If is then natural to define 
the capacity of a node of degree k as 

qk = (w)k+ma k , (2) 

where to is a real positive parameter and characterizes 
the tolerance of the node to load. A node is considered 
to fail if the number of walkers w on the node exceeds this 
capacity. Therefore, the probability F\y 0 (k) that a node 
of degree k experiences an overload failure is calculated 
by summing up the distribution function hk{w) over w 
larger than qk . Using the regularized incomplete beta 
function I x (a,b) for this summation [62j, the overload 
probability is expressed as [58 ] 

F\y 0 (k) = h/2M 0 i\qk\ + 1 , W 0 — ( 3 ) 

where the floor function [x\ represents the greatest inte¬ 
ger less than or equal to x. 

Applying the above idea of the overload probability, a 
cascade of failures starting with a specific large network 
can be described as follows [57j . 

(i) Prepare an initial connected, undirected, and un¬ 
correlated network Q 0 with N 0 nodes and M 0 edges, 
in which Wq random walkers exist, where Wq is 
chosen so as to be proportional to M 0 . In addition, 
determine the capacity qk of each node by Eq. ©. 

(ii) At each time step r, assign W T random walkers to 
the network Q T at step r, where the total load W T 
is given by 

w ' = w °{w) T - < 4 > 

Here M T is the total number of edges in the network 
Q t and r is a real positive parameter. 

(iii) Calculate the overload probability of every node, 
and remove nodes from Q T with this probability. 

(iv) Repeat (ii) and (iii) until no node is removed in the 
procedure (iii). 

The reduction of the total load in the procedure (ii) cor¬ 
responds to actual cascades of failures during which the 
total load is reduced to some extent to prevent the break¬ 
down of the network function. We call the exponent r in 
Eq. Jp the load reduction parameter hereafter. It should 
be emphasized that the overload probability in the pro¬ 
cedure (iii) cannot be calculated by Eq. (0 with W 0 re¬ 
placed by W T for the following two reasons. First, the 
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degree fc of a node in the network Q T is not the same with 
its initial degree k 0 . The capacity of a node is determined 
by the initial degree fco, while the probability to find a 
walker on this node is proportional to the present de¬ 
gree k. Thus, the overload probability depends on both 
k and ko. Secondly, the network Q T is not necessarily 
connected, though the initial network Go is connected. If 
Q T is not connected, W T random walkers are distributed 
to each component in proportion to the number of edges 
in the component, before starting the next cascade step. 
Taking into account these remarks and the fact that ran¬ 
dom walkers cannot jump to other components, the over¬ 
load probability of a node of degree fc, whose initial de¬ 
gree is fco, in the a-th component of Q T is given by 

F w? ( k 0 , k) = I k/ 2 M? (L<?fc 0 J + 1, W“ - Lo'fcoJ), (5) 

where is the number of edges in the a-th component 
of Q t and W“ = W T M™/M T is the load assigned to the 
a-th component. Since k is always equal to fco and the 
network Q T is connected at r = 0, Fw?{ko,k) at r = 0 
coincides with Fw 0 (k) given by Eq. 0 . Thus, Eq. 0 is 
a general form of the overload probability for r > 0. 

The relative size Sf of the giant component in the net¬ 
work Qf at the final stage of the cascade process is an 
important quantity to evaluate the robustness of net¬ 
works against cascading overload failures. This quan¬ 
tity Sf can be calculated by combining the generating 
function method icfl and the master equation for the 
probability n r (fco,fc) that a node in Q T has the present 
degree fc and the initial degree fco, without simulating 
numerically the cascade process (i)-(iv) [57j . By means 
of this method, it has been clarified that there exists a 
threshold value of the load reduction parameter r c (7V 0 ) 
above which Sf becomes finite and below which Sf = 0 
and r c (N(j) for No —> oo provides a percolation transition 
by cascading overload failures [5^]. The critical property 
of Qf at r = r c has also been confirmed by the fractality 
of the giant component in Gf ■ These facts will be closely 
related to the results of the present work. However, in 
the case that the structure and node capacities of a net¬ 
work with which the cascade starts depend on results of 
cascading failures occurring in the past, as in the model 
of this work, it is unfortunately impossible to apply the 
method utilizing the generating function. In such a case, 
we need to simulate numerically the process (i)-(iv) to 
find the final network state after a cascade. 


B. Network evolution 

The basic idea of our dynamical model is to combine 
the growth of a network and its decay into smaller ones 
by cascading overload failures. The overload probability 
F\v 0 (fc) given by Eq. 0 is almost independent of M 0 
if Wq oc Mo- However, the probability that the first 
failures inducing a cascade of subsequent failures occur 
in the network increases with the network growth. Thus, 
we expect that the network cannot grow infinitely and its 


size fluctuates around a certain value. Since our model 
includes two different types of dynamics, i.e., growth and 
cascade of failures, in order to distinguish them clearly, a 
time step of the network growth is hereafter denoted by 
t in parentheses and a cascade step by subscript r. The 
concrete algorithm of the network evolution in our model 
is then given as follows: 

(1) Start with a small and connected network G{ 0) with 

iVini nodes and edges, in which Wi n i random 
walkers exist. Wj n i is set as W- ln i = where a 

is a positive constant. The capacities of the nodes 
in £7(0) are calculated by Eq. (0. 

(2) At every time step t > 1, add a new node with /x 
edges, where /x is in the range of 2 < /x < N n i, and 
connect the new node to /x different nodes selected 
randomly from the network G(t — 1) at time t — 1. 
Let Go(t) be the network at this stage. 

(3) Place Wo(t) = aMo{t ) random walkers on the net¬ 
work Go(t), where Mo(t) is the number of edges in 
Go(t), and calculate the capacity of the new node 
by using Eq. 0 with fc replaced by /x and for the 
total load Wo(t). 

(4) Perform cascading overload failures starting with 
Go {t) in accordance with the process (i)-(iv) de¬ 
scribed in Sect. Ill Al In this cascade process, iso¬ 
lated zero-degree nodes generated by the elimina¬ 
tion of all their adjacent nodes, in addition to the 
overloaded nodes themselves, are removed from the 
system. Let G(t) be the resultant network after 
completing the cascade. 

(5) Repeat the procedure from (2) to (4) for a suffi¬ 
ciently long period. 

We should make several remarks concerning the above 
algorithm. In the procedure (2), the number of edges ^ 
of a newly added node must be larger than 2. Other¬ 
wise, once the network is divided into disconnected com¬ 
ponents by cascading failures, components never merge. 
The network after a long time becomes an assembly of a 
large number of small graphs. Thus, the dynamical sys¬ 
tem has a qualitatively different property from that for 
H > 2. It should be also emphasized that the network 
Go(t) is not necessarily connected. If Go(t) consists of 
plural components, the total load Wo(t) is distributed to 
each component in the same way as the case of cascad¬ 
ing failures. Namely, the number of walkers in the a-th 
component is allocated by 


W 0 a (t) 


W 0 (t) 


mi m 

Mo (t ) 


aM 0 “(f), 


( 6 ) 


where is the number of edges in the a-th compo¬ 

nent of Go (t). The calculation of the new-node capacity 
in the procedure (3) is actually done by using Wo(t), be¬ 
cause random walkers cannot move to other components. 
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The capacity of the node that is introduced at time t and 
embedded in the a-th component is then presented by 

<f(t) = (™>M + ( 7 ) 


where ( w) fl = Wg{t)p“ = a/i/2, = yj(w)^{1 - p^) , 

and = fj,/2Mg(t). Once the node capacity is deter¬ 
mined, the value of q never changes until the node is 
eliminated. 

The process of cascading overload failures in the pro¬ 
cedure (4) basically follows the steps (i)-(iv) in Sect. Ill Al 
with replacing Go, N 0 , M 0 , and W 0 by Go(t), N 0 (t), 
and W 0 (t), respectively, where N 0 (t) is the num¬ 
ber of nodes in Go(t). In addition to the possibility of 
Go(t) being disconnected, there are two other differences 
in the detailed treatment of the cascade process. Firstly, 
the overload probability of a node cannot be written as 
Eq. @. This is because the node capacity given by 
Eq. 0 depends on when the node was introduced in 
the system. Thus, the overload probability of the node i 
at cascade step r is written as 

F w?(t)(i) = (IVOOJ + 1 ,W?(t) - \_q P (U)\) , 

( 8 ) 

where fc,; is the degree of the node i, ti is the time at 
which the node i was introduced, a and /? are the indices 
of the components to which the node i belongs at the 
present cascade step r and at ti, respectively, and q^(ti) 
is presented by Eq. 0. The symbols W“(t) and M“(i) 
represent the number of random walkers and number of 
edges in the a-th component of G r (t), where G T (t) is the 
network at the step r of the cascade starting with Go ( t ). 
The second difference is in the load reduction scheme dur¬ 
ing the cascade. In Sect. Ill Al the total number of random 
walkers is reduced in accordance with the reduction of the 
network size during the cascade, as expressed by Eq. 0, 
to prevent the breakdown of the network function. It is 
actually difficult to reduce quickly the total load when the 
network size becomes large. Taking into account such re¬ 
alistic situations of cascading failures, the load reduction 
parameter r characterizing how quickly the total load is 
reduced with the reduction of the network size should 
decrease with No(t). Therefore, the total load during the 
cascade is reduced according to 


W T (t) = W 0 (t) 


M r jt) 

M 0 (t) 


r l N o (*)] 


(9) 


where M T (t) is the number of edges in Grit) and r[iVo(t)] 
is a decreasing function of No(t). Since large-scale cas¬ 
cades are more likely to occur if r is small, the property 
that r[iVo(t)] decreases with No(t) prevents the network 
from growing too large. 

If all nodes are eliminated from the system in the pro¬ 
cedure (4), the network G(t) is reset to £(0) and con¬ 
tinue the network evolution from the procedure (2). In 
the present model, the time scale of cascades measured 
by the step r is assumed to be much faster than that of 


the network growth measured by the step t, and the re¬ 
laxation time of random walkers in a network is further 
shorter than a single cascade step. We concentrate, in 
this work, on the temporal evolution of networks in the 
time scale of the network growth. Therefore, information 
on the network G(t) at every end of the procedure (4) is 
recorded to investigate the model. 


III. RESULTS AND DISCUSSION 


Our model includes several parameters and conditions. 
These are the numbers of nodes N- m \ and edges M\ n \ in 
the initial network 1/(0), the topology of 17(0), the load 
carried by a single edge a, the node tolerance parameter 
to, the number of edges of a newly added node p, and 
the functional form of the load reduction parameter r(N). 
In this section, we fix these parameters as follows: The 
initial network 17(0) is a triangular ring with A r ; n i = 3 
and Mi n i = 3. The parameters a and /i are set as a = 2.0 
and /i = 2, respectively. The value of to is chosen from 
the range of 5.0 < m < 7.0. The function r(N) is set as 


r(N) 


^max 

T max (M max - N) 


M max - N-, 


ini 


for 2 < N < iVjni, 
for Mini <N< M max , 
for N > M max . 

( 10 ) 


This function decreases from its maximum value r max to 
zero as N increases. Since a cascade of overload failures 
with r = 0 eliminates all nodes in any network, M max 
gives a rough estimation of the maximum network size 
in the dynamics. Here we set M max = 1, 000 and r max = 
1.0. We will explain later the reason why we adopt the 
above parameter values and discuss suitable ranges for 
the model parameters to obtain SOC character in the 
network evolution. 


A. Number of nodes and other network measures 

We first examine the time dependence of the size of 
the network G(t). Figure 1 shows the number of nodes 
N(t) in G(t) for the first 10 4 time steps. The result clearly 
demonstrates that the network cannot grow infinitely and 
the size N(t) largely fluctuates by repetitive growth and 
decay of the network. We can find some features in the 
line shape of N{t). In the early stage, N(t) increases 
almost monotonically with time, because the probability 
that any of the nodes in the network fail is low due to 
small No(t) [= N(t — 1) + 1]. Namely, in this time region, 
the expectation number of failed nodes is less than 1. Af¬ 
ter this region, the expectation value becomes larger than 
1 and some nodes fail in Go(t). However, since r[M 0 (t)] for 
still small No(t) is rather large, cascades are not widely 
spread. Thus, N(t) for 400 < t < 800 keeps increasing 
with relatively small drops. When N(t) becomes larger 
than 800, r[IVo(i)] is so small that a cascade of failures 
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FIG. 1. (Color online) Time dependences of the number of 
nodes N(t) in the network Q(t) (thick black line) and number 
of nodes in its largest component (thin red line). The 

node tolerance parameter is set as m = 5.0. 
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never stops until all nodes are eliminated. Such complete 
collapses occur at t = 1, 824, 4,603, and 6, 374 in Fig. |T| 
After a complete collapse, the system evolves in a simi¬ 
lar manner to the evolution from t = 0. In addition to 
N(t), we plot in Fig. [T] the size of the largest component 
Ni,c(t) contained in Q{t). The size IVlcW basically fol¬ 
lows the variation of N(t) at most of the time steps, but 
sometimes iV L c(t) drops substantially though N(t) does 
not change so much. At these times, the network is de¬ 
composed into small components by cascading overload 
failures. The statistics of magnitudes of drops in N(t), 
i.e., cascade sizes, will be argued in the next subsection. 

We also calculated several quantities that character¬ 
ize the network structure at time t. Figure [2] shows the 
average degree (fc), the clustering coefficient C, and the 
average path length (l) of the network Q(t) as a func¬ 
tion of t. These quantities for the largest component of 
Q(t) take almost the same values as those for Q(t). The 
average degree fluctuates around (k) = 2(= p) though 
it becomes significantly larger than this value immedi¬ 
ately after a complete collapse. We have confirmed that 
the degree distribution V{k) of Q(t) hardly depends on 
time and decays exponentially for large k if N(t) is large 
enough (not shown here). This is reasonable because 
random attachment of new nodes and cascading failures 
without introducing degree correlations make the net¬ 
work topology similar to a homogeneous random graph. 
The clustering coefficient C is quite small (C )$ 10 -4 ), 
except for Q(i) at and just after complete collapses. At 
a complete collapse, C is equal to 1 because Q(t) is a 
triangular ring at this time. The clustering coefficient, 
however, rapidly decreases with the network growth by 
random attachments. Sometimes C becomes equal to 
zero, which implies that the network takes a tree (or for- 


FIG. 2. Time dependences of the average degree (top), the 
clustering coefficient (middle), and the average path length 
(bottom) of the network G(t). The node tolerance parameter 
is set as m = 5.0. These quantities for the largest component 
of Q(t) are not shown in this figure, because their line shapes 
almost overlap with those for G(t). 

est) structure. The average path length ( l } of the net¬ 
work at a complete collapse is obviously 1, and after that 
( l ) increases gradually with relatively large fluctuations. 
Considering that N(t) is less than 1,000, ( l ) close to or 
more than 10 is too large to regard the network as being 
small world. Then, we can expect that Q{t) giving very 
large ( l ) has a fractal structure. Before discussing the 
fractality of generated networks, it will be examined in 
the next subsection whether our dynamical system ex¬ 
hibits SOC behavior. 


B. Avalanche size and self-organized criticality 

Sudden drops of the network size N(t) found in Fig. [T] 
corresponds to decays of the network by cascading over¬ 
load failures. Magnitudes of these drops represent scales 
of cascading overload failures. Here we define the 
avalanche size S(t) as the number of nodes that are re¬ 
moved during a single cascade of overload failures occur¬ 
ring at the time t. Figure [3] shows the avalanche size 
S(t) obtained from N(t) shown in Fig.Q] The avalanche 
size largely fluctuates even if one ignores huge S(ty s at 
complete collapses. Values of S(t) at most of the time 
steps are less than 50, while on rare occasions S(t) ex¬ 
ceeds 300. The inset of Fig. [3] demonstrates that these 
avalanches occur intermittently with inactive intervals. 
This intermittency suggests a possibility that the net- 
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FIG. 3. Time dependence of the avalanche size S(t) obtained 
from N(t) shown in Fig. [T] To make the figure easy to see, 
cascades leading to complete collapses are eliminated from 
the figure. The intermittency of cascades is found in the inset 
that magnifies the main figure for 100 <t< 300. 



FIG. 4. Distribution function W ( T ) of the inactive time in¬ 
terval T during which no overload failure occurs. The distri¬ 
bution W ( T ) is obtained from the dynamics up to t = 5 x 10 6 
under the condition m = 5.0. The inset shows the distribution 
of N(t) in the same dynamics. 


work dynamics possesses SOC characteristics. In order 
to find further evidences of SOC dynamics, we exam¬ 
ine the distribution function W(T) of the inactive time 
interval T between avalanches. The distribution W(T) 
obtained from the dynamics under the same conditions as 
those for Fig. Q] but continued up to 5 x 10 6 time steps is 
presented in Fig. Q] This figure clearly shows that W(T) 




FIG. 5. (Color online) (a) Distribution function P(S) of the 
avalanche size S and (b) the distribution function n(s) of 
the cluster size s for different values of the node tolerance 
parameter, i.e., m = 5.0, 6.0, and 7.0. The distributions P(S) 
and n(s) are obtained from the dynamics up to t = 5 x 10 6 . 
In both panels, the distributions for m = 6.0 and 7.0 are 
vertically shifted for clarity. 


obeys a power law, 

W(T) oc T _T? , (11) 

in an intermediate region of T. The least-squares fit for 
the data within 20 < T <50 gives ry = 3.00 ± 0.03. The 
small hump near T ~ 150 comes from a finite-size effect 
related to the existence of the most probable network 
size ATtyp. This size is about 700 for our choice of the 
model parameters as depicted in the inset of Fig. 2] We 
have confirmed the correlations between N{t) and S(t) 
and between the inactive interval T after a cascade and 
the cascade (avalanche) size S. These correlations and 
iVtyp ~ 700 lead the frequently-appearing time interval 
at T ~ 150. 

The distributions of the avalanche size S for several 
values of the node tolerance parameter m are presented 
in Fig. (5ja) . The avalanche size distribution P(S) also 
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follows a power-law relation, i.e., 

P{S) oc S~ x . (12) 

Our result indicates that the exponent A does not depend 
on to and is estimated as A = 2.60 ± 0.02 for m = 7.0. 
The rightmost hump in P(S) at around S ~ 900 repre¬ 
sents the contribution from complete collapses. The sec¬ 
ond hump from the right corresponds to decays of large 
networks to assemblies of dimers. Moreover, the broad 
hump near S = 150 found in the result for m = 5.0 is re¬ 
lated to the hump in W(T) shown by Fig. [4] This broad 
hump represents the typical avalanche size of cascades 
from typical networks with N typ ~ 700 nodes. There¬ 
fore, these humps are attributable to finite-size effects 
associated with our choice of the model parameters. 

Figure (3Jb) shows the cluster size distribution n(s) for 
three different values of to. The cluster size s at time t 
is the number of nodes in a component included in the 
network Q{t). The cluster size distribution function n(s) 
is calculated from all components of the network at every 
time step in the entire dynamics. The distribution n{s ) 
has a power-law form, 


n(s) oc s T , (13) 

as well as W(T) and P(S). The exponent r, calculated 
as r = 2.92 ± 0.02 for to = 7.0, is also independent of 
to. In contrast to the distribution P(S), the influence of 
complete collapses to n(s) is inconspicuous. This is be¬ 
cause the size of a network just after a complete collapse 
is ./Vin;(= 3) and the number of components with s = 3 
generated in the whole period of the network evolution 
is extremely large compared to the number of complete 
collapses occurring in the same period. 

All the above results, namely the intermittency of 
S(t) and the power-law forms of W(T), P{S), and n(s), 
strongly support that the dynamics of network structure 
in our model exhibits SOC behavior. These results also 
show that the universality class of self-organized critical¬ 
ity does not depend on the node tolerance parameter to. 
The relation to other parameters will be discussed later. 


C. Fractal and small-world networks 

As we mentioned in Sect. nm a cascade of overload 
failures gives a fatal damage to a connected network 
Go of size A 0 if the load reduction parameter r is less 
than r c (No), and the giant component after a cascade at 
r = r c (A 0 ) has a fractal structure .57'|. In the present 
SOC model, on the other hand, r decreases with the net¬ 
work size Na(t). For r[A 0 (f)] chosen as Eq. (flOl) . the 
parameter r decreases from a large enough value r max 
for No(t) < A; n ; to zero for A 0 (t) > A max . Since any 
network is completely collapsed by a cascade of failures 
at r = 0, through the network growth, r[Ao(f)] must 
eventually encounter the critical value r c at which the 
cascade of failures provides the percolation transition of 



Network size 


FIG. 6. Histogram of the number of pre-critical networks 
just before critical cascades versus the size of the pre-critical 
network. The histogram is obtained for the first 10 4 critical 
cascades occurring in the dynamics under the condition m = 
5.0. 


the network. (Precisely speaking, the term “critical” is 
not appropriate because No(t) is finite. However, we use 
this terminology by supposing the case that sufficiently 
large networks are generated under suitable values of the 
model parameters. The word “critical” in the rest of 
this paper will be used in the same sense.) If a network 
with the size N 0 (t) satisfying r[A 0 (f)] = r c experiences 
a cascade of overload failures, we expect that the giant 
component after the cascade has a fractal structure. 

In the case of cascading failures starting with a fixed 
connected network Go of size Aa in which all the capac¬ 
ities of nodes are definitely determined by their degrees 
and the initial total load Wo, the value of r c (Ao) is theo¬ 
retically calculated as addressed in the last paragraph of 
Sect. Ill Al 57(. In our SOC model, however, the capacity 
of a node depends on the total load at the time when 
the node was introduced in the system. Thus, the node 
capacities in the network Go(t) depends strongly on the 
past history of Qo(t), and the critical load reduction pa¬ 
rameter r c cannot be uniquely determined by the size of 
Go{t)- Since the theoretical method proposed by Ref. ITVl 
is not applicable to dynamics governed by such hysteresis 
effects, we need to examine numerically whether r[Aa(t)] 
of the network Go(t ) is close to an unknown value of r c 
peculiar to Go(t). 

If r[Ao(t)] of the network Go{t) is much larger than r c , 
a cascade of overload failures, if any, eliminates only a 
small fraction of nodes from Go(t) and does not change 
the giant component size so much. On the other hand, a 
cascade with r[Ao(f)] <C r c causes a complete collapse of 
Go{t)- If r[Ao(i)] is close to r c , the cascade is marginal, 
for which the size of the giant component after the cas¬ 
cade must be much smaller than the original giant com¬ 
ponent size of Go(t) but still much larger than A; n ;. From 
the above consideration, we regard in this work a cascade 
of overload failures at time t satisfying the following con- 
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FIG. 7. (Color online) Nb(Ib) for giant components in crit¬ 
ical networks generated by SOC dynamics under the condi¬ 
tions m = 5.0, 6.0, and 7.0. The longitudinal axis indicates 
Nb(Ib)/Nb(1) averaged over 1, 000 realizations of critical net¬ 
works. The results for m = 6.0 and 7.0 are vertically shifted 
for clarity. The straight dashed line has the slope ds = 1-53 
which is obtained by the least-squares fit for the data for 
m = 7.0 from Zb = 2 to 20. 

ditions as a critical cascade whose load reduction param¬ 
eter should be close to r c : 

Nictt-l) ~ 0,5 ^ NhC ^ ~ 100, (14) 

where -/Vlc(Z) is the number of nodes in the largest com¬ 
ponent oiQ(t). The specific values 0.5 and 100 in Eq. (fldl) 
are not important as long as Alc (t)/NhC (t — 1) and 
TVlc (t) are much smaller and larger than 1, respectively. 
In the sense of the percolation transition by cascading 
overload failures, a network after completing a critical 
cascade can be considered as a critical network Q c . Also, 
we call a network just before a critical cascade a pre- 
critical network C/ pre . Figure [6] shows the histogram of 
the number of pre-critical networks as a function of the 
size of f/pre- This result indicates that critical cascades 
are likely to occur on networks of size N(t) ~ 700 for the 
present parameter set. Considering that the most prob¬ 
able network size iV typ is also about 700 as shown by 
the inset of Fig. [2 critical cascades take place frequently 
during SOC dynamics. This means that critical networks 
are generated very often by such cascades. 

We study the fractal property of giant components 
in critical networks. As we explained in Sect. [I] if a 
given connected network is fractal, the minimum number 
Nb(Ib) of subgraphs of diameter less than Ib required to 
cover the network satisfies the relation 

JVb(Zb) (15) 

where dp is the fractal dimension of the network Q. 
We calculate JVb(Zb) for giant components in 1,000 


FIG. 8. (Color online) Nb(Ib) for largest components in net¬ 
works that first reach the size N(t) = 300 after complete col¬ 
lapses in SOC dynamics under the conditions m = 5.0, 6.0, 
and 7.0. The longitudinal axis indicates Nb(Ib)/Nb(1) aver¬ 
aged over 1, 000 realizations of such networks. The results for 
m = 6.0 and 7.0 are vertically shifted for clarity. The dashed 
lines are guides to the eye. 


critical networks appearing in the dynamics by using 
the compact-box-burning algorithm [22| and average 
-^b(Zb)/ZVb( 1) over th ese realizations, where .ZVb(I) is 
equal to the number of nodes in the largest component. 
The results for m = 5.0, 6.0, and 7.0 are plotted in Fig. [7] 
These plots clearly demonstrate that the quantity IVb(Zb) 
satisfies Eq. m and the fractal dimension ds does not 
depend on m. The value of de estimated from the re¬ 
sult for m = 7.0 is 1.53 ± 0.01. It is interesting that 
this fractal dimension is close to ds = 1.54 ± 0.01 that 
has been computed for the giant component after a crit¬ 
ical cascade starting with an Erdos-Renyi (ER) random 
graph (57). The topology of a pre-critical network Q pre 
in SOC dynamics is not the same as that of the ER ran¬ 
dom graph 1 /er- In addition, the capacity of a node in 
t/pre depends on the total load of the system when the 
node was introduced, while the node capacity in C/er is 
determined by the degree of the node and a fixed ini¬ 
tial total load Wq. In spite of these discrepancies, it is 
not surprising that both fractal dimensions are the same, 
which implies the same universality class between perco¬ 
lation transitions for Q pie and C/er- This is because the 
degree distribution of Q pre has an exponential tail, like 
that of <7er, and the node capacity distribution is not 
wide, which behaves similarly to the distribution of N(t) 
shown in the inset of Fig. [IJ 

Let us examine TVb(Zb) for off-critical networks. If a 
network is far from criticality, we can expect that the 
network has a small-world structure, because the net¬ 
work formed by random attachment of new nodes has 
many short-cut edges. For a small-world network, the 
number of covering subgraphs ZVb(Zr) decreases expo- 





9 


‘4 



FIG. 9. (Color online) (a) Time dependence of the number of 
nodes Ni,c(t) in the largest component of the network Q{t) in 
the dynamics for m = 7.0. The arrows indicate the times at 
which N b (IbYs are calculated, (b) Ab(Zb)’s for largest com¬ 
ponents in Q(t i), G(t 2 ), G(tz), and G{tY) from top to bottom. 
Thin lines are guides to the eye. The rapid decrease of JVb(Zb) 
for Ib 3> lco at t = ts or t± indicates that the largest compo¬ 
nent has a small-world structure in a longer length-scale than 
lco- The inset shows ZVb(Zb) at t = t 4 in a semi-logarithmic 
scale. 


nentially with Zb, namely, 

N b (Ib) oc exp(-Z B /Z 0 ), (16) 

where Iq is a characteristic path length. We calculated 
ZVb(Zb) for networks (or their largest components if not 
connected) that first reach the size N(t) = 300 after com¬ 
plete collapses, and averaged over 1,000 

realizations of such networks in SOC dynamics. The re¬ 
sults shown in Fig. [ 8 ] indicate the small-world property 
of these networks. Our SOC model thus generates both 
fractal and small-world networks in a single dynamics. 

We further investigate the crossover behavior from 
fractal to small-world structure associated with the time 
evolution from a critical network. Figure [9] illustrates a 
typical profile change of Nb{Ib) for networks formed at 
several times from t\ at which a critical network appears 


to ti just before the next complete collapse. The times at 
which Nb(Ib)'s are calculated are indicated in Fig. H a ) 
by arrows on the time dependence of the largest com¬ 
ponent size iVLc(i)- At t = t\, N b (Ib) follows a power 
law, which suggests that the giant component in G(t\) 
has a fractal structure as we expect. After this time, 
the largest component size rapidly increases as shown in 
Fig. Ha). This is because newly added nodes are more 
likely to merge separated fractal components but less 
likely to be connected onto a single component. There¬ 
fore, the largest component at t = O remains fractal at 
almost any scale. When the time elapses further, the in¬ 
crease of N^cit) becomes moderate. This implies that 
the merging process of separated components has been 
mostly finished and new nodes are simply incorporated in 
the largest component. In this case, newly added nodes 
bring short-cut edges in the largest component, which 
makes the network small-world as shown by Nb{1 B ) at 
t = £3 in Fig.Hb). More precisely, the network is small- 
world in a longer length scale than the average distance 
l co between terminal nodes of short-cut edges introduced 
by new nodes, while it is fractal for Zb -C l c o- This sit¬ 
uation is similar to the case that a lattice-like network 
changes into a small-world one by random rewirings in 
the Watts-Strogatz model [64] In fact, a high density of 
short-cut edges at t = reduces the crossover length l co 
and the small-world property can be found in the whole 
Zb range as shown by the inset of Fig. Hb). 


D. Suitable choice of parameter values 

All the above arguments are based on specific values 
of the model parameters. If the dynamical properties 
presented above are peculiar to these parameter values, 
it cannot be said that the present model exhibits self- 
organized criticality, because of the necessity of tuning 
the external parameters. It has, however, been confirmed 
that the results are essentially independent of the choice 
of parameter values if these parameters lie in suitable 
ranges. In this subsection, we discuss the suitable pa¬ 
rameter ranges to realize SOC dynamics. 

To find the suitable ranges of parameter values, let us 
consider how large a network could grow if the system 
did not experience any critical cascades and complete 
collapses. Even in this case, a network cannot grow in¬ 
finitely. The expectation number ( S ) of eliminated nodes 
per unit time step increases with the network size ZV(Z), 
and eventually ( S ) reaches the incrementation of N ( t ) at 
every time step due to the participation of a new node. 
Once this is the case, the network does not grow any 
more. The network size N(t ) then fluctuates around a 
stationary size with small amplitudes. The stationary 
size N st can be roughly estimated by the overload prob¬ 
ability. In the absence of critical cascades and complete 
collapses, we can consider approximately that all cascad¬ 
ing overload failures stop at the first step of the cascade 
process and subsequent avalanches triggered by the first 
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failures do not occur, because avalanche sizes are small. 
This approximation enables us to calculate the steady- 
state expectation number of failed nodes per unit time 
step by 


(S) = N st J2Vst(k)F W[Net) (k), (17) 

k 

where V s t (k) is the degree distribution of a steady-state 
network Q st and F W (n s t )(k) is the overload probability of 
a node of degree k in Q st . With the aid of the regularized 
incomplete beta function, the probability F w ^ Bt )(k) is, 
with reference to Eq. ©, given by 


F W (N st) (k) = I k/ltNmt [[q^N st )\ + 1, W(iV st ) - [q^N st )\], 

(18) 

where 


q^(N s t) — 


afj, 



and 


W(N st ) 


afj,N st 

2 


(19) 


( 20 ) 


Here, we approximated the average degree of Q s t by /i 
for the reason that nodes with degree greater than [i are 
more likely to be eliminated by overload failures in C/ St 
while the average degree would be larger than // (equal 
to 2/i) if the network monotonically grew without any 
node elimination. In the steady state, ( S) must be equal 
to the incrementation of the network size per unit time 
step, namely 1. Therefore, the stationary size N st is de¬ 
termined by the relation, 


N at 


1 

J2k 'Pst(k)F W ( Nst) (k) 


( 21 ) 


It is thus significant to elucidate how N st and N pre de¬ 
pend on a (the load carried by a single edge), in (the 
node tolerance parameter), /i (the number of edges of a 
newly added node), and the functional form of the load 
reduction parameter r(N). 

The stationary size N st depends on a, m, and /. i . Equa¬ 
tion m reveals the relation of N s t to these parameters. 
Since the preferential elimination of nodes with degree 
much larger than // in Gst gives a sharp peak of V s t (fc) at 
k = fi, V s t {k) hardly depends on N s t . Furthermore, the 
overload probability Fyy^^^k) presented by Eq. (1151) in¬ 
deed depends only very weakly on iV st , which comes from 
the property of the regularized incomplete beta func¬ 
tion. Therefore, Eq. m is not actually transcendental, 
and N st can be evaluated by V s t{k) and F W ( N ){k) for 
a haphazardly chosen value of N{~^> /i). The probability 
F W (N){k) given by Eqs. (fl8l) - (l20l) with N st replaced by N 
is a decreasing function of to and // for any k. Hence, N st 
obtained by Eq. m increases with m and /i regardless 
of Vstik). The a dependence of EV(Af)(fc) is, however, 
influenced by the value of k. Fyy^^k) increases with a 
if k > fx, while it decreases for k < fi. Meanwhile, for a 
fixed value of a, F^^){k) for k < // is negligibly small. 
Thus, -FV(jv) ( k ) for k larger than p dominates the sum¬ 
mation in Eq. ED, independently of the form of Vst(k). 
This fact and the property of Fyy^^k) of being an in¬ 
creasing function of a for k > // show that N s t decreases 
with a. Consequently, we need to choose large values of 
to and n and a small value of a to obtain large N s t- 

On the other hand, N pie is the typical size of a network 
whose load reduction parameter r is equal to the criti¬ 
cal value r c specific to the network. The parameter r is 
uniquely determined by the network size N, while r c de¬ 
pends not only on N but also on the past history of the 
network. Approximating the typical size of pre-critical 
networks by the size of a typical pre-critical network Q pre , 
N pie must satisfy 


If we neglect critical cascades and complete collapses, the 
network can grow up to the size N st obtained by solving 
the above transcendental equation. But actually, the net¬ 
work encounters critical cascades or complete collapses 
before reaching this size if N s t is larger than the typical 
size N pie of pre-critical networks C/ pre . In this case, and 
only in this case, critical cascades generate fractal net¬ 
works, and the present model exhibits SOC character. 
Otherwise, critical cascades themselves never take place 
in the dynamics. Thus, the condition to realize SOC dy¬ 
namics is 

Nst > N pre > 1, (22) 

where the condition N pie 1 guarantees that the system 
is large enough to exhibit genuine self-organized critical¬ 
ity. What is the relation between the above condition 
and the model parameters? Among several parameters 
characterizing our model, parameters related to the ini¬ 
tial network G(0), namely, A r j n i, M; n i, and the topology 
of 1 /( 0 ), are obviously irrelevant to the condition (E^l) . 


r(N pie ) = r c [G pie (N pie )} : (23) 

where r c [t/pre(A r p re )] is the critical load reduction param¬ 
eter of C/pre whose size is N pre . Since the right-hand size 
of Eq. (1^51) is a function of a, to, /k, and Ap re , the size 
N pie as the solution of Eq. (l23ll depends on these param¬ 
eters in addition to the functional form of r(N). If r(N) 
decreases slowly with N, however, the solution A^ pre is 
mainly governed by the form of r(N) rather than the 
precise value of r c . In order to satisfy A^ pre 1, r(N) 
needs to decrease very slowly with the network size. In 
the case that r(N) is set as Eq. (fTUl) . N max must be cho¬ 
sen to be large enough. 

In conclusion, the load reduction parameter r(N) must 
decrease with N very slowly to realize the condition 
Npi e ^ 1 , and the node tolerance parameter m and the 
load by edge a should be large and small enough, re¬ 
spectively, so that N s t becomes much larger than N P i e . 
Although a large value of /i is preferable for the condi¬ 
tion N s t N pre , results are not strongly influenced by 
/i because the number of edges of a new node is always 
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restricted by 2 < /i < A^ n ; with small N - ln Our choice of 
values for a, m, and /i in this section obviously satisfies 
the condition N st N pre , because we have critical cas¬ 
cades in the dynamics. In fact, N st estimated by Eq. m 
with the numerically obtained V s t(k) is 861 for a = 2 . 0 , 
/r = 2, and m = 5.0, which is larger than N pie ~ 700 as 
indicated in Fig [ 6 ] We have confirmed that the univer¬ 
sality class of self-organized criticality, namely the set of 
the exponents 77 , A, r, and g?b, does not depend on the 
choice of parameter values if the condition (B^l) is satis¬ 
fied. It has also been checked that the functional form of 
r(N) is irrelevant to SOC dynamics as far as r(N) is a 
slowly decreasing function of N. 

IV. SUMMARY 

We have proposed a model of self-organized critical 
(SOC) dynamics of complex networks and presented a 
possible explanation of the emergence of fractal and 
small-world networks. Our model combines a network 
growth and its decay due to the instability of large grown 
networks against cascading overload failures. Cascad¬ 
ing failures occur intermittently and prevent networks 
from growing infinitely. The distribution of the inactive 
time interval between successive cascades of failures has a 
power-law form. Both the avalanche size that is the num¬ 
ber of eliminated nodes in a single cascade and the cluster 
size defined as the number of nodes in a connected com¬ 
ponent also obey power-law distributions. These facts 
indicate that the network dynamics possesses SOC char¬ 
acteristics. During the SOC dynamics, the load reduc¬ 
tion parameter r varies with the network size. When 
r of the network coincides with its critical value r c , a 
cascade of overload failures (critical cascade) decays the 
network into a critical one. We have shown that giant 
components just after critical cascades have fractal struc¬ 
tures. The fractal dimension de is close to that for the 
giant component after a critical cascade starting with 
an Erdos-Renyi random graph. In contrast, networks 


far from criticality display the small-world property. In 
particular, we demonstrated the crossover behavior from 
fractal to small-world structure in a growing process from 
a critical network, which is caused by short-cut edges in¬ 
troduced by newly added nodes. We have also discussed 
suitable parameter values to realize SOC dynamics. 

It is significant to notice that the present model is 
somewhat different from previous SOC models. In a 
conventional SOC model, a routine procedure in the dy¬ 
namics, such as placement of grains of sand in the sand- 
pile model [25j or renewals of fitness values in the Bak- 
Sneppen model J2fj, takes a system close to the criti¬ 
cal point, but the instability of critical or near-critical 
states drives the system away from criticality accompa¬ 
nied by some sort of avalanches. In our model, on the 
other hand, the network growth as a routine procedure 
takes the system away from the critical point, but the 
instability of large grown networks makes the network 
critical. Although the roles of growth and instability are 
opposite to those of conventional models, the system de¬ 
scribed by our model exhibits the most of SOC charac¬ 
teristics as explained in Sect. m This implies that the 
present model provides a new type of self-organized criti¬ 
cality. Our model generates non-scale-free networks with 
homogeneous degree distributions and belongs to a spe¬ 
cific universality class of SOC dynamics, independently 
of the choice of values of the model parameters. It is then 
interesting to study how the model should be modified to 
belong to another class of self-organized criticality with 
forming scale-free networks. 
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